%     VIRTUAL LAMINATE - программа для получения эффективных свойств слоистых пластин
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES, ELEMENTS, RVE_Coordinates, Node_Merge_Tolerance, Ply_Elements] = Generate_Mesh(Materials, Stacking, Element_Order, N_Elements_Ply, N_Elements_Side, Element_Critical_Ratio)

global Merge_Tolerance_Multiplier;

fprintf (1, 'Построение сетки...                                     ');

% Минимальная толщина слоя
t_min = Materials{Stacking(1,3)}.t;

for i=1:size(Stacking,1)
    
    if (Materials{Stacking(i,3)}.t < t_min)
        t_min = Materials{Stacking(i,3)}.t;
    end
    
end

% Минимальная высота элемента
h_min = t_min / N_Elements_Ply;

% Длина и ширина элемента
Element_Length = h_min * Element_Critical_Ratio;

% Длина и ширина ячейки периодичности
RVE_Length = Element_Length * N_Elements_Side;

% Построение плоской сетки
[NODES_2D, ELEMENTS_2D, Next_Node_Number, Next_Element_Number] = Mesh_2D_Quadrilateral([0,0], [RVE_Length,0], [0,RVE_Length], [RVE_Length,RVE_Length], N_Elements_Side, N_Elements_Side, 1, 1, 0); %#ok<ASGLU>

% Длина минимального ребра элемента
Min_Element_Edge = h_min;

% Координаты узлов по высоте
Node_Heights = zeros(10000,1);
Node_Heights(1,1) = 0;
counter = 1;

Ply_Numbers = zeros(10000,1);
counter_2 = 0;

Ply_Elements = zeros(size(Stacking,1),1);

for i=1:size(Stacking,1)
    
    t_ply = Materials{Stacking(i,3)}.t;
    
    h_elm = t_ply / N_Elements_Ply;
    
    Base_Height = Node_Heights(counter,1);
    
    Element_Ratio = h_elm / Element_Length;
    
    if (Element_Ratio > Element_Critical_Ratio)
        
        N_Elm_local = ceil(N_Elements_Ply * Element_Ratio / Element_Critical_Ratio);
        
        for j=1:N_Elm_local
            
            counter = counter + 1;
            Node_Heights(counter,1) = Base_Height + t_ply * j / N_Elm_local;
            
            counter_2 = counter_2 + 1;
            Ply_Numbers(counter_2) = i;
            
        end
        
        Ply_Elements(i) = N_Elements_Side * N_Elements_Side * N_Elm_local;
        
    else
        
        for j=1:N_Elements_Ply
            
            counter = counter + 1;
            Node_Heights(counter,1) = Base_Height + t_ply * j / N_Elements_Ply;
            
            counter_2 = counter_2 + 1;
            Ply_Numbers(counter_2) = i;
            
        end
        
        Ply_Elements(i) = N_Elements_Side * N_Elements_Side * N_Elements_Ply;
        
    end
    
end

Node_Heights = Node_Heights(1:counter,:);

Ply_Numbers = Ply_Numbers(1:counter_2,:);

% Число элементов по высоте ламината
N_Elements_High = size(Node_Heights,1) - 1;

% Допуск при слиянии узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Число узлов плоской сетки
N_Nodes_2D = size(NODES_2D,1);

% Выделение памяти под узлы объемной сетки
NODES = zeros((N_Elements_High + 1) * N_Nodes_2D, 4);

% Заполнение массива узлов объемной сетки
for i=0:N_Elements_High
    
    NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D,:) = NODES_2D(:,:);
    
    if (i > 0)
        
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 1) = NODES((i-1)*N_Nodes_2D + 1:i*N_Nodes_2D, 1) + ones(N_Nodes_2D,1) * N_Nodes_2D;
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 4) = ones(N_Nodes_2D,1) * Node_Heights(i+1);
        
    end
    
end

% Число элементов плоской сетки
N_Elements_2D = size(ELEMENTS_2D,1);

% Выделение памяти элементы объемной сетки
ELEMENTS = zeros(N_Elements_High * N_Elements_2D, 12);

% Заполнение первого слоя объемной сетки
for i=1:N_Elements_2D
    
    switch ELEMENTS_2D(i,2)
        
        case 1
            
            ELEMENTS(i,1) = ELEMENTS_2D(i,1);
            ELEMENTS(i,2) = 2;
            ELEMENTS(i,3) = 4;
            ELEMENTS(i,4) = ELEMENTS_2D(i,4);
            
            ELEMENTS(i,5) = ELEMENTS_2D(i,6);
            ELEMENTS(i,6) = ELEMENTS_2D(i,5);
            ELEMENTS(i,7) = ELEMENTS_2D(i,5) + N_Nodes_2D;
            ELEMENTS(i,8) = ELEMENTS_2D(i,6) + N_Nodes_2D;
            
        case 2
            
            switch ELEMENTS_2D(i,3)
                
                case 3
                    
                    ELEMENTS(i,1) = ELEMENTS_2D(i,1);
                    ELEMENTS(i,2) = 3;
                    ELEMENTS(i,3) = 6;
                    ELEMENTS(i,4) = ELEMENTS_2D(i,4);
                    
                    Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,7),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));
                    
                    if (dot(Element_Normal, [0,0,1]) > 0)
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        
                    else
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        
                    end
                    
                case 4
                    
                    ELEMENTS(i,1) = ELEMENTS_2D(i,1);
                    ELEMENTS(i,2) = 3;
                    ELEMENTS(i,3) = 8;
                    ELEMENTS(i,4) = ELEMENTS_2D(i,4);
                    
                    Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,8),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));

                    if (dot(Element_Normal, [0,0,1]) > 0)
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,8);
                        ELEMENTS(i,9) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,11) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,12) = ELEMENTS_2D(i,8) + N_Nodes_2D;
                        
                    else
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,8);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,9) = ELEMENTS_2D(i,8) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,11) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,12) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        
                    end
                
            end
        
    end
    
end

% Заполнение остальных слоев объемной сетки
Mask = ELEMENTS(1 : N_Elements_2D, 5:12) > 0;
Mask = double(Mask);
    
for i=1:N_Elements_High-1
    
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 1) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 1) + ones(N_Elements_2D,1)*N_Elements_2D;
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 2:4) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 2:4);
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 5:12) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 5:12) + Mask*N_Nodes_2D;
    
end

% Номера слоев
for i=1:N_Elements_High
    
    ELEMENTS(1 + (i-1)*N_Elements_2D : i*N_Elements_2D, 4) = ones(N_Elements_2D,1) * Ply_Numbers(i);
    
end

% Габариты ячейки
RVE_Coordinates = zeros(3,2);

RVE_Coordinates(1,1) = 0;
RVE_Coordinates(1,2) = RVE_Length;
RVE_Coordinates(2,1) = 0;
RVE_Coordinates(2,2) = RVE_Length;
RVE_Coordinates(3,1) = 0;
RVE_Coordinates(3,2) = Node_Heights(size(Node_Heights,1));

fprintf (1, 'OK!\n');

% Конвертация элементов в элементы второго порядка
if (strcmp(Element_Order, 'QUADRATIC'))
    
    fprintf (1, 'Конвертация элементов в квадратичные...                 ');
    
    % Уменьшение допуска для слияния узлов
    Node_Merge_Tolerance = Node_Merge_Tolerance/2;
    
    [NODES, ELEMENTS] = Convert_Mesh_To_Quad(NODES, ELEMENTS, Node_Merge_Tolerance);
    
    fprintf (1, 'OK!\n');
    
end
